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Hydrodynamic turbulence driven by crust-core differential rotation imposes a fundamental noise 
floor on gravitational wave observations of neutron stars. The gravitational wave emission peaks 
at the Kolmogorov decoherence frequency which, for reasonable values of the crust-core shear, Af2, 
occurs near the most sensitive part of the frequency band for ground-based, long-baseline interfer- 
ometers. We calculate the energy density spectrum of the stochastic gravitational wave background 
from a cosmological population of turbulent neutron stars generalising previous calculations for in- 
dividual sources. The spectrum resembles a piecewise power law, Q sv ,(v) = Q a v a , with a = — 1 
and 7 above and below the decoherence frequency respectively, and its normalisation scales as 
Q. a oc (AO,) 7 . Non-detection of a stochastic signal by Initial LIGO implies an upper limit on Af2 
and hence by implication on the internal relaxation time-scale for the crust and core to come into 
co-rotation, r c i — Afi/f2, where Cl is the observed electromagnetic spin-down rate, with Td < 10 7 yr 
for accreting millisecond pulsars and Td < 10 J yr for radio-loud pulsars. Target limits on Td are also 
estimated for future detectors, namely Advanced LIGO and the Einstein Telescope, and are found 
to be astrophysically interesting. 

PACS numbers: 95.85.Sz 04.30.Db 97.60.Jd 



I. INTRODUCTION 



The electromagnetic braking torque acting on a neu- 
tron star drives persistent differential rotation between 
the rigid crust and the multiple superfluid components 
in the interior, causing observed phenomena like ro- 
tational glitches [H-dj]. Classical relaxation processes, 
like Ekman pumping and Sweet-Parker circulation, and 
quantum mechanical relaxation processes, like superfluid 
vortex creep, determine the long-term angular velocity 
shear between the various internal components [343- I n 
turn, differential rotation drives turbulence when the 
Reynolds number is high (Re > 10 11 ) as in a neutron 
star [1,0]. The turbulence takes two distinct forms: (1) a 
Kolmogorov-like cascade of macroscopic "eddies" or cir- 
culation cells 043 i an d (2) a self-sustaining tangle of 
microscopic quantised vortices flH [l6j , created when the 
rectilinear vortex array in a uniformly rotating superfluid 
is disrupted by instabilities driven by meridional circula- 
tio n Ha . [3 , interfacial and bulk two-stream instabilities 
and nuclear pinning forces [20L HH . 



Turbulence driven by differential rotation is ax- 
isymmctric when averaged over long times but non- 
axisymmetric instantaneously It therefore emits stochas- 
tic gravitational radiation. The root-mean-square gravi- 
tational wave strain assuming incompressibility and long- 



term isotropy is given by Q 
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where M+ and i?* are respectively the mass and radius 
of the neutron star, d is the distance from Earth and 
Ail is the angular velocity difference between the slower 
crust (angular velocity il) and faster core. The emission 
is predicted to be astrophysically relevant [and poten- 
tially detectable by third-generation interferometers like 
the Einstein Telescope (ET)] for sources including pro- 
toneutron stars [Hj], accreting millisecond pulsars, ac- 
creting white dwarfs on the verge of accretion-induced 
collapse [3 and young pulsars with super-rotating 
cores, whose deceleration is inhibited by buoyancy Q. 
Even before any gravitational-wave detections, the indi- 
rect spin-down limit from radio timing observations puts 
interesting upper limits on Ail (see the left panel of fig- 
ure 5 in Ref. @), with Ail/il < 10~ 2 in some sources, 
approaching the shear inferred from glitch data. 

In this paper we extend the single-source calculations 
in Ref. to calculate the stochastic gravitational wave 
signal coming from all neutron stars in the Universe. We 
then use existing non-detections by the Laser Interferom- 
eter Gravitational- wave Observatory (LIGO) to place an 
upper limit on Ail across the cosmological population as 
a whole, generalising the single-source limits. Many kinds 
of sources contribute to the stochastic background in 
the LIGO frequency band; for reviews of various sources 
see Refs. j25M27| . Searches from the first generation of 
LIGO-class detectors have already pushed the gravita- 
tional wave energy density, il gw (u), at frequency v, below 
that inferred from Big Bang nucleosynthesis (BBN) and 
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the cosmic microwave background |28[ . Cross-correlation 
searches are performed assuming a power law of the form 
= fl a v a ; they place an upper bound on Q a given 
q [29 - l32j . We show in this paper that the gravitational- 
wave energy density from turbulent neutron stars is well 
approximated by a piecewise power law with a rising com- 
ponent (a = 7) at low v and a decaying tail (a = —1) 
at high v. The power laws join and peak near the tur- 
bulence decoherence frequency, with f2 gw (j/) at the peak 
scaling with the seventh power of Aft. 

The structure of the paper is as follows. In section [TT1 
we review the relevant statistical properties of supcrfluid 
turbulence and calculate the gravitational wave signal 
from a single turbulent neutron star pi BP and multiple 
stars pi C[) . In section UTTl we use LIGO non-detections 
to derive upper limits on Af2 under three scenarios: a 
universal A 17 set by nuclear physics pil Al) . a broad A17 
distribution set by the balance between internal damp- 
ing and electromagnetic spin down B ip , and a narrow 
A17 distribution for the population of accreting millisec- 
ond pulsars pil B 2p . The prospects for direct detection 
of this cosmological background by current and future 
detectors are summarised in section IPV1 



II. GRAVITATIONAL RADIATION FROM 
NEUTRON STAR TURBULENCE 

A. Turbulence statistics 

We first review briefly the main statistical properties 
(e.g. the autocorrelation function) of the stochastic grav- 
itational wave signal emitted by supcrfluid turbulence in 
a single, differentially rotating neutron star. We refer 
the reader to Ref. for details of the derivation. The 
stochastic background resulting from the superposition 
of multiple sources at multiple redshifts is calculated in 
section III CI 

The wave strain in the transverse-traceless gauge as 
measured by an observer at a distance d is given by 
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where S is the (£, m) current multipole written as a 
function of the retarded time, t, and T-? 2 ' m is the tensor 
spherical harmonic describing the angular dependence of 
the radiation field (see equation 2.30f of Ref. |33j). In the 
Newtonian approximation (i.e. slow internal motions and 
weak internal gravity), which is adequate for describing 
subsonic turbulence driven by slow differential rotation 
(i.e. A17/17 <C 1), the current multipole is 
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where Y lm is a scalar spherical harmonic, v(x, t) is the 
turbulent velocity field and p(x) is the fluid mass den- 
sity. For incompressible turbulence, the mass multipolcs 
vanish to a good approximation, as the turbulent mo- 
tions are subsonic and p(x) is uniform 1 . The gravita- 
tional wave signal is dominated by I = 2 for most real- 
istic neutron stars, i.e. objects with Reynolds number 
Re< (c/R+An) 8 [|. 

We consider stationary, isotropic turbulence, for which 
the mean wave strain at the observer vanishes, and the 
leading non-zero moment is the autocorrelation function, 
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with r = t' — t. where (. . .) denotes the ensemble average 
over realisations of the turbulence. From equation ([3]), 
the mean-square wave strain, h^. ms = C(0), is propor- 
tional to the zero-lag autocorrelation of the fluctuating 
vorticity field, curlv(x, t). 

Three-dimensional, global simulations of shear-driven 
neutron star turbulence suggest that, for a two- 
component Hall-Vinen-Bekarevich-Khalatnikov supcr- 
fluid, the flow is approximately isotropic and station- 
ary for Re > 10 4 0, [IMIj]. We adopt the standard 
Kraichnan-Kolmogorov form for the unequal-time veloc- 
ity correlator [1, [Hi 111 2 

(v m (k, t) v q (k', if)*) = VP mq (k)F (k, t) 6 (k - k') . (5) 

Here, V is the total volume of the system, P mg (fc) = 
&mq — k m k q is a, projection operator perpendicular to the 
wave vector k, P(k) is the power spectral density [i.e. 
(27r)~ 3 fc 2 P(fc) is the kinetic energy in the flow per unit 
wavenumber] , and we have 



F(k,r) = P(yfc)exp [-7ny(fc) 2 T 2 /4] 
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where r](k)~ l is the eddy turnover time at wavenumber 
fc, and e is the energy dissipation rate per unit enthalpy. 

The power spectral density is not known for supcrfluid 
turbulence. In the absence of specific information, we as- 
sume the Kolmogorov law, P(k) oc k a , with a = —11/3, 
which characterises isotropic, high-Re, Navier-Stokes tur- 
bulence. The power law extends between the stirring 
wavenumber, k s = 27r/i?*, where i?* is the stellar ra- 
dius, and the viscous dissipation wavenumber, kd = 



1 Density perturbations are of order pR^ (Af2) 2 /c 2 , where c s is the 
sound speed. Hence the mass quadrupolc is smaller by a factor 
RirAQc/Cg than the current quadrupole. See \Si for recent high- 
resolution simulations of compressible, relativistic turbulence. 

2 Equation J5)l differs from equation (2) of Ref. @| by a factor of 
(2tt) 3 , which we have absorbed into the definition of the power 
spectral density, P(k). 
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[8e/ (27w 3 )] 1/4 , where v is the kinematic viscosity. The 
Kolmogorov law is consistent with preliminary numerical 
simulations in the context of the two-component Hall- 
Vinen-Bekarevich-Khalatnikov model with 10 4 < Re < 
10 5 in the viscous component 0. However, P(k) remains 
unknown in both stratified Navier-Stokes turbulence and 
unstratified supcrfluid turbulence even in terrestrial ex- 
periments [33| • The subtle spatio-temporal anisotropies 
caused by layering, intermittency and rotational polari- 
sation (e.g. of the superfluid vortex tangle) are discussed 
briefly in Appendix [SJ with pointers to the voluminous 
literature discussing these issues. 

Combining the above equations, and expanding the 
plane-wave Fourier components (|5|) in spherical harmon- 
ics to evaluate S 2m , one arrives at the following formula 
for the autocorrelation function for the mode (£, m) = 
(2,m): 
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In Ref. [| only the I = m = 2 mode was considered. 
In reality all the modes I = 2, |m| < 2 contribute to the 
wave strain. For a single source, the angular dependence 
of the combined wave strain is complicated. However one 
can show that all the modes £ = 2, |m| < 2 contribute 
equally to the total autocorrelation function, which is 
the quantity of interest for a stochastic background from 
many isotropically distributed sources; see equation 
below. 

The decoherence time corresponding to the half-strain 
point, C(t c ) = h 2 . ms /2, is 



r c = 0.35r ? (fc s )- 1 
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where Af2 = e 1 / 3 i?* 2 ^ 3 is the angular velocity lag be- 
tween crust and core. 

Equation is derived assuming \T^ k 2 ' 2m \ 2 = 1, which 
is true only for an optimally situated observer; that is, 
C(r) in equation ([8]) is the maximum of equation (j4]) over 
observer orientation for a given m. As we are ultimately 
interested in the .stochastic background from an isotropic 
population of cosmological sources (section III C|) . it is 
more appropriate to replace \T^ k 2 ' 2m \ 2 = 1 with the sky- 
averaged tensor harmonic product 
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summed over the repeated indices j, k, where <p and 9 
represent latitude and longitude of the Earth relative to 
the source. Hence equation ([8|) must be multiplied by 
5/(47r) in what follows. 



B. Single sources 

A key ingredient in calculating the stochastic back- 
ground from multiple sources is the energy spectrum 
emitted by a single source. To calculate the radiated 
energy per unit area per unit frequency in terms of C (r) , 
we begin with the standard form of the radiated energy 
per unit area per unit time expressed in the transverse- 
traceless gauge [e.g. HU 
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The ensemble average in equation f(12|) is over realisations 
of the turbulence, as in equation Q. Integrating both 
sides of equation ([12"]) with respect to retarded time, we 
apply Parseval's theorem to obtain [39| 
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where T [. . .] denotes the Fourier transform. The above 
quantity is evaluated in the emitted (i.e. comoving) 
frame, remembering that the sources are at rcdshift 
z > 0. To clearly distinguish the comoving and observers 
frames, we denote the frequency emitted (observed) with 
(without) a subscript 'e', such that v e ~ v{\ + z). 

The Wiener-Khintchine theorem relates the spectral 
density to the inverse Fourier transform of the autocor- 
relation function (eg., [!(|)- Employing the identity (4l| 
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valid for any stationary, differcntiable random variable, 
X(t), we obtain 
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where T is the emitting lifetime. The emitted energy in- 
creases monotonically with time while the source is on. 
In equation (|15|) , the Wiener-Khintchine theorem is eval- 
uated in the limit where T is finite but much greater than 
r c ; the details are presented in Appendix [Bj 

Twice differentiating the autocorrelation function 
given by equation ((5]) and taking the Fourier transform 
according to equation (|15[) . we reach the final result 




emitted during this time interval equals 
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FIG. 1: Gravitational wave energy per unit area per unit 
frequency emitted by a single neutron star as a function of 
frequency, for various values of the angular velocity shear AS1. 



with rjs — f]{k s ). This function is plotted in normalised 
form in figure Q] for a typical neutron star with A/* = 
1.4 Mq, i?* = 10 km, kinematic viscosity, v = lm 2 s _1 
and for various Af2 values 3 . The spectrum peaks near 
the inverse of the decoherence time, r c , given by equation 
(fTTJ)) . From equations ([1]), (fTU)) and (|TT5| . one can show 
that the peak energy per unit area per unit frequency 



scales as Ail 7 , with d 2 E K 
and oc v~ 2 for v e > l/r c . 
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C. Multiple sources 

To calculate the total energy emitted by multiple, con- 
tinuously emitting sources at multiple redshifts z e , con- 
sider an infinitesimal time interval dt(z e ) = (dt / dz e )dz e 
between lookback times t(z + dz e ) and t(z e ). The en- 
ergy per unit area per unit logarithmic frequency interval 



3 We scale v to 1 m 2 s _1 , instead of the standard 10m 2 s — 1 [^be- 
cause Landau damping by weakly screened transverse plasmons 
[43l lowers v by nearly an order of magnitude relative to the 
standard electron-electron scattering value. The viscosity only 
affects the highest and lowest frequencies, and therefore plays an 
insignificant role in determining the integrated stochastic back- 
ground. 
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where d 2 E gw /dSdi/ e is given by equation (|15|). with T 
replaced by dt(z e ), the factor (1 + z e ) _1 accounts for the 
redshiftcd energy of the gravitons, and v e = v(\ + z e ). 
The number density of sources emitting during the time 
interval equals the total number of neutron stars born at 
redshifts z b > z e , 
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where p^z^) is the star formation rate per unit comoving 
volume, $(Af) is the initial mass function, (M m j n , M max ) 
defines the initial mass function range that forms neutron 
stars, M* = 1.4 M© is the neutron star mass, and dt/dz b 
is set by the cosmology. Throughout the article we adopt 
a concordance cosmology with Sl m = 0.26, f^A = 0.74 
and H = 73kms _1 Mpc" 1 . We allow for the possibility 
that d 2 E gw /dSdv e is also a function of z b through Afl, 
which decreases as the source spins down and is there- 
fore a function of its age. Following convention, we can 
then express the energy density in the gravitational wave 
stochastic background, fraction of the closure 

energy density per logarithmic frequency interval, 
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(20) 



In flU), p c c 2 = 3H$c?/8nG is the critical energy density 
required to close the Universe and v is the frequency in 
the observer's frame. 

Superficially, equation (fhSj) and (|20|) look identical to 
equation (5) of Ref. (lij . Physically, however, the fac- 
tors in the integrand have different interpretations, be- 
cause our sources emit continuously, while those in Ref. 
phi are discrete, short-lived bursts (i.e. inspirals). In 
Ref. [Hj, N(z e )dz e equals the infinitesimal number of 
inspiral events occurring between z e and z e + dz e , while 
dEgw/dQxii/e) is the total, non-infinitesimal energy per 
logarithmic frequency emitted by each event. In con- 
trast, in (|20[) . N(z e ) is the non-infinitesimal number of 
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continuously emitting neutron stars in existence between 
z e and z e + dz e , while dE gw / d(\nv e ) equals the infinitesi- 
mal energy per logarithmic frequency emitted during the 
time interval dt(z e ) = (dt/ dz e )dz e . 

From equations (|16p and (f2"0")l , one can show that the 
stochastic background is described approximately by a 
piecewise power law, Q gvf {v) = Q a v a , with a = 7 and 
a = — 1 for v < v c and v > v c respectively. The peak 
frequency, f c , is a population- weighted average of the re- 
ciprocal of the decoherence time given in (fTU)) (see section 



For the remainder of the article we adopt the modified 
Salpeter A initial mass function and the corresponding 
parametric fit for the star formation rate from Hopkins 
and Beacom (4f|. For safety, we also verify our calcu- 
lations against the star formation rate given in Cucciati 
et al. |4g|; the results are similar, with Q gw (v c ) differ- 
ing by < 6% between the two mass functions. As a 
rule, we take the range of zero-age main sequence masses 
that form neutron stars to be 8 < Mj M@ < 40. Vary- 
ing the minimum and maximum masses over the ranges 
4 < M min /M© < 8 and 20 < M max /M < 40 respec- 
tively changes Q gw {v c ) by < 25%, one of the smaller un- 
certainties in our overall calculation. The initial mass 
range for neutron star formation is discussed more fully 
inRef. El 



III. MAXIMUM SHEAR FROM 
GRAVITATIONAL WAVE NON-DETECTIONS 

The peak gravitational wave energy density, Q g -w(v c )i 
from individual neutron stars scales as the seventh power 
of AQ as shown in section IHBI However, detailed first- 
principles predictions for AQ are not yet available, nor 
is there any compelling observational support (e.g. from 
radio pulsar timing) for any particular choice of AQ. For 
the moment, therefore, we are obliged to use gravita- 
tional wave non-detections to place an upper limit on 
AQ across the neutron star population. We do this for 
three illustrative astrophysical scenarios in this section: 

• a unique AQ value across the population, as ex- 
pected if the shear is set by the balance between 
Magnus and nuclear pinning forces in the inner 
crust superfluid; 

• a broad AQ distribution, wherein the shear is 
approximately proportional to the observed spin- 
down rate, balanced by some sort of internal relax- 
ation (e.g. vortex creep, viscous damping or Sweet- 
Parker circulation); and 

• a narrow AQ distribution at AQ > 10 rad s -1 , due 
to fast accretion-driven spin-up in accreting mil- 
lisecond rotators like low-mass X-ray binaries. 
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FIG. 2: Gravitational wave energy density, Sl gw , as a function 
of frequency, u, for three values of the angular velocity shear, 
AQ. = lOrads -1 (solid blue curve), 30rads _1 (dotted green 
curve) and 50rads _1 (dashed red curve). The three solid 
(dashed) black curves are the one- (three-) year noise curves 
for Initial LIGO, Advanced LIGO and ET, running from top 
to bottom. 



A. Unique shear 

We begin by calculating flg W (v) assuming a constant 
shear, AO, in every neutron star in the Universe. In fig- 
ure [H we plot f2 gw as a function of v for Afl = 10, 30 
and 50 rad s~ . We also show the gravitational wave de- 
tection limits for Initial LIGO, Advanced LIGO and the 
proposed Einstein Telescope (ET) , each for one and three 
years of data collection (thick solid and dashed black 
curves res pec tively). The sensitivity curves are taken 
from Ref. [48| , where an analytic fitting formula is given 
for the noise power spectral density of each detector (see 
table 1 of Ref. HH). The conversion to il gw for two co- 
located detectors and uncorrelated instrumental noise is 
given in section 8.1.2 of Ref. (48|. Throughout this arti- 
cle we define a signal as being detectable if the predicted 
signal amplitude lies above the noise curve for that par- 
ticular instrument at any frequency. In reality, a more 
rigorous cross-correlation analysis will need to be done 
to detect a source. 



For 10 rad s -1 < AQ < 10 2 rads _1 , fl gw (v) peaks near 
the most sensitive part of the LIGO frequency band. The 
strong scaling of O gw (^ c ) cx AO 7 is clear in figure [5] and 
implies that the stochastic background from neutron star 
turbulence may be observable by ET for AQ in the above 
range. From current non-detections with Initial LIGO, 
we obtain AQ < 55rads _1 . A hypothetical Advanced 
LIGO non-detection with three years of data would imply 
AQ < 25rads _1 , which is high but not unphysically so. 
The latter limit is competitive with the indirect spin- 
down limit AQ/Q < 0.04 inferred from radio timing Q 
for that subset of the neutron star population with Q > 
7.5 x 10 2 rads _1 . 
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B. Shear distributions 

Let us now assume that the shear is proportional to 
the spin-down rate, where the constant of proportionality 
equals the internal relaxation time-scale, Td- We consider 
two neutron star populations, one based on the observed 
spin-down rates of radio-loud pulsars, the other based on 
the observed spin-up rates of X-ray accreting millisecond 
pulsars. 



1. Radio-loud pulsars 

From the distribution of radio-loud pulsars in the 
ATNF catalogue [4£| 4 , we select objects with surface 
magnetic field greater than 2 x 10 10 G to exclude mil- 
lisecond pulsars, which are treated in section UlI B 21 We 
fit a log-normal distribution 5 to Aft = T^Cl, leaving Td as 
a model parameter to be constrained by future (current) 
gravitational wave (non)-detections. Wc assume Td is the 
same in all objects for simplicity and is set by the viscos- 
ity for example. Note that the Afl distribution must fall 
off faster than AO 7 otherwise fl gw diverges. Hence we 
cut off the log-normal distributions at a maximum shear 
Af2 max , related to the centrifugal break-up angular ve- 
locity, a conservative choice. 

In figure [3] we plot ft gw as a function of v for Td = 10 8 
and 10 10 s. The solid, dotted and dashed curves cor- 
respond to Afi max /27r = 2, 1.5, 1 kHz respectively. A 
comparison of figures [2] and [3] shows that the background 
from the radio-loud distribution is dominated by objects 
with AO near Af2 max ; as Af2 max decreases, Q gw (v) turns 
over at lower frequencies and fl gv/ (u c ) decreases. 

The total gravitational wave energy density, ob- 
tained by integrating the curve in figure |3J fl g ^ = 
Jq 00 dlni/Og W (^), places additional constraints on 
^gw(^)- For example, fl^ must be smaller than the 
baryon energy density inferred from cosmological obser- 
vations, i.e. < fib ~ 0.04 [Hi]]; clearly the back- 
ground studied in this paper, which is emitted primarily 
at z ~ 1, ultimately comes from mechanical energy in 
baryons created at higher z. This leads to upper limits 
on Td that are comparable with those derived in figure 
[3J for the ET sensitivity curve with Afi max /27r ss 2 kHz 
and for Advanced LIGO with ASl max /27r ps 1 kHz. 



4 http: / / www.atnf.csiro.au/peoplc / pulsar/ psreat / 

5 The decision to fit a log-normal distribution to the data was 
motivated by inspection (given other uncertainties, qualitative 
agreement is sufficient), as well as by population synthesis mod- 
els that fit, for example, log-normals to the magnetic field (B) 
distribution |5C|. which is related to Q through B 2 Q S oc Q. We 
quantify how good this fit is to the data by calculating the ratio 
of the moments of the data to that of the functional fit. We find 
that the ratios corresponding to the first and second moments 
are both unity by construction, while the third and fourth are 
respectively 0.979 and 0.94f . 
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FIG. 3: Gravitational wave energy density assuming the Afl 
distribution for all neutron stars is identical to the Afl dis- 
tribution for radio-loud pulsars, with Ail = TdQ, where Td 
is an internal relaxation time-scale and the Cl distribution 
is drawn from the ATNF catalogue [4jJ. The solid, dotted 
and dashed coloured curves correspond to maximum shear 
Afi ma x/27r = 2, 1.5 and 1kHz respectively. The black curves 
are the sensitivity curves for Initial LIGO, Advanced LIGO 
and ET assuming observation times of one year (solid curves) 
and three years (dashed curves). 



We emphasise however, that large shears, i.e. 
AS1/27T > 1 kHz, are unlikely astrophysically, except 
perhaps in a small subset of young objects with super- 
rotating cores Q. Hence, the lesson of figure[3Jis not that 
detections are expected at v > 1 kHz from a handful of 
strongly sheared pulsars, but rather that (i) upper lim- 
its on Afl are best obtained from 0.2 kHz < v < 1 kHz, 
where the theoretical curves hug the detector sensitivity 
curves; and (ii) such upper limits on Afl are approxi- 
mately independent of Af2 max . 

Figure [3J is drawn assuming that all neutron stars 
are described by the radio-loud distribution. In real- 
ity, only a certain fraction, Af are described by this 
distribution; most neutron stars lie beyond the pul- 
sar death line, where magnetospheric electron-positron 
pair cascades switch off [e.g. [52j. From (|2T))) . we have 
Qgw (A/ - , v) = Affl gw (1, v), and the curves in figure [3J are 
upper limits; in general they are lower by the factor Af. 
In figure|3]we plot the Initial LIGO, Advanced LIGO and 
ET non-detection curves as a function of Td and Af. That 
is, for any ordered pair (r<j, Af) that lies above the curve 
plotted in figure 21 some portion of f2 gw (^) lies above 
the sensitivity curve of that particular detector configu- 
ration. The shaded region above the blue line in figure 
S] indicates the parameter space that has already been 
ruled out by the non-detection of a stochastic signal by 
Initial LIGO. 

For any given value of Td, Advanced LIGO provides a 
limit on Af two orders of magnitude better than Initial 
LIGO. The improvement with ET is three orders of mag- 
nitude. Alternatively, for realistic neutron star popula- 
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FIG. 4: Non-detection limits from Initial LIGO, Advanced 
LIGO and ET in the (jd, VV)-plane, where Td is the internal 
relaxation time-scale and AT is the fraction of neutron stars 
following the AS! distribution for radio-loud pulsars. A one- 
year observation is assumed. The absence of detection by 
LIGO so far rules out the parameter space above the blue 
line. 

tions with 10~ 3 < M < 1(T 5 , Initial LIGO non-detection 
implies Td < 10 13 s. The limit becomes Td < 5 x 10 11 s 
and Td < 10 10 s for Advanced LIGO and ET respectively. 
We comment briefly on how close these limits are to chal- 
lenging theoretical expectations in section HVl 

One can set a complementary upper limit on the com- 
bination of Td and Af from the condition 51*^: < f2f, dis- 
cussed above. Again, for A£7 max /27r w 2kHz (1 kHz), 
this upper limit is similar to the ET (Advanced LIGO) 
curve presented in figure |4] We reiterate that such large 
shears are astrophysically very unlikely. Future genera- 
tions of gravitational wave detectors with higher sensi- 
tivity will place more interesting limits. 



2. Accreting millisecond pulsars 

The scaling Jl gw (^ c ) oc (AJ1) 7 means that relatively 
few sources with large shears can dominate the back- 
ground. One natural place to find large shears are ac- 
creting millisecond pulsars. The shear is maintained by 
the accretion spin- up torque, A^ acc , which is often greater 
than an isolated star's electromagnetic spin-down torque 
[Hi . Typically we have N acc = 111 m M y/GM+R*, where 
M is the mass accretion rate, / is the moment of inertia 
and the lever arm (i.e. magnetospheric radius) is approx- 
imately R+. There are only a few X-ray timing observa- 
tions of f2 [e.g., see[13], and we therefore retain AJ1 = r f /fi 
as the free variable that we wish to constrain. Similar to 
section Till B 1| Jl gw is a function of two parameters: the 
millisecond pulsar fraction A/", and r^fi. In figure [S] we 
plot the non-detection curves for the Initial LIGO, Ad- 
vanced LIGO and ET detectors in the (t^O, A/)-plane. 



FIG. 5: Non-detection limits from Initial LIGO, Advanced 
LIGO and ET for the population of accreting millisecond sys- 
tems in the (rdfi, A/")-plane, where AQ — rati is set by the 
accretion spin-up torque (see text) and Af is the fraction of 
neutron stars that are accreting millisecond pulsars. A one- 
year observation time is assumed. The absence of detection 
by LIGO so far rules out the parameter space above the blue 
line. 

Non-detection by Initial LIGO implies the shaded blue 
region is already ruled out. 

From Ref. |53 |. a typical value of M for accret- 
ing millisecond pulsars leads to Cl ~ 10~ 12 s~ 2 . Pop- 
ulation synthesis models suggest there are ~ 10 3 low- 
and intermediate-mass X-ray binaries out of ~ 10 9 com- 
pact objects in the galaxy (e.g., [H, [5(| and references 
therein) . Taking this as typical of the universal popula- 
tion, i.e. Af ~ 10~ 6 , Initial LIGO non-detection implies 
T d < 10 15 s. Advanced LIGO and ET will push this limit 
to Td < 10 14 s and Td < 10 13 s respectively. 

IV. CONCLUSION 

In this paper we have calculated the stochastic grav- 
itational wave background from supcrfluid turbulence 
driven by differential rotation in a cosmological popula- 
tions of neutron stars, generalising the single-source cal- 
culation in Rcf. We found that the gravitational 
wave energy density per logarithmic frequency interval 
peaks at v c ~ AJ1 and that its peak value scales as 
Sl gw (f c ) oc (AJl) 7 . Hence, relatively few sources with 
large shears dominates the background. 

We evaluated the background in three specific sce- 
narios. Firstly, we took all sources to have a unique 
A£7, as when balance between the Magnus and nu- 
clear pinning forces sets the shear in the inner crust 
superfluid. It was found that the background is de- 
tectable by third-generation ground-based gravitational 
wave detectors such as the proposed Einstein Telescope 
for AO > 20 rad s" 1 . Secondly, we took the known dis- 
tribution of radio-loud pulsars from the ATNF catalogue 
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|49( to be representative of all neutron stars, the shear 
was then assumed to be proportional to the spin-down 
rate, where the constant of proportionality is the relax- 
ation time-scale for the core and the crust to come into 
co-rotation in the absence of spin-down. These models 
are then parametrised by the fraction of total objects that 
lie within this distribution, Af. For Af ~ 10~ 3 , current 
LIGO non-detection implies < 1 x 10 13 s. Advanced 
LIGO and ET non-detection limits are tj<5x 10 11 s 
and 1 x 10 10 s respectively. Finally, we calculated the 
background from accreting, rapidly rotating systems like 
low-mass X-ray binaries. Assuming again a common re- 
laxation time-scale, f2 = 2n x 10~ n rads~ 2 from X-ray 
timing data, and a reasonable value of Af ~ 10~ 6 , non- 
detection by LIGO implies Td < 1 x 10 15 s. Advanced 
LIGO and ET will push the limit to r d < 1 x 10 14 s and 
1 x 10 13 s respectively. 

The gravitational wave energy density approximately 
follows a piecewise power law, Q gw = Sl a v a , with a = 7 
for v < v c and a = — 1 for v > v c . LIGO and Virgo 
cross-correlation searches for a stochastic background 
have looked for power laws, albeit with — 3 < a < 3 
[28l [29l [3lL HH ; see also Ref. 57} for a discussion of pa- 
rameter estimation within these models. It will be worth 
extending the range of exponents to include a = 7 in the 
future. 

There are a number of uncertainties in our results that 
require further investigation. In Ref. @, it was shown 
that the root-mean-square wave strain, h Tms , is insen- 
sitive to the exponent of the turbulence power spectral 
density. On the other hand, the form of the eddy turnover 
time, 7y(A:) _1 , and hence the form of the autocorrelation 
function, C(r), depends more closely on the dynamics of 
the turbulent eddies, and these dynamics are poorly un- 
derstood in superfluid turbulence [37||- Superfluid spher- 
ical Couette simulations modelling neutron star turbu- 
lence exhibit a Kolmogorov-likc cascade 

[SEMI, but 

they are filtered spectrally to ensure numerical stability, 
so more work needs to be done before their output can 
be trusted fully. In addition, stratification is likely to 
modify the turbulent spectrum (see Appendix [XJ, and 
magnetic fields must eventually be incorporated too (see 
discussion in Ref. Q). 

In light of the above uncertainties, the results in this 
paper are deliberately expressed as upper limits on Af2 
from non-detections, rather than predictions on fl gw (v) 
given a known AQ. It is of course pertinent to ask how 
the upper limits on Afl compare with the best guesses 
for Af2 in the literature. One can approach this question 
in many ways. 

• In glitches, the observed fractional angular velocity 
jump is 1CT 11 < An/n < 10~ 4 , although the ab- 
sence of a reservoir effect, whereby the glitch size is 
proportional to the time elapsed since the preced- 
ing glitch, suggests this change in angular velocity 
is a small fraction of the underlying shear 0, HH . 

• If the angular velocity lag between the crust and 



the superfluid core is set by balance between the 
Magnus and pinning forces, the differential ang ular 
velocity can be as large as Af2 ~ lrads -1 [21j . 

• Gravitational wave emission from hydrodynamic 
turbulence removes rotational kinetic energy from 
a neutron star, causing the star to spin down. A 
fundamental upper limit follows from noting that 
this gravitational wave spin-down is less than the 
spin-down observed in radio timing experiments, 
giving AQ/Q < 10" 2 [|. 

• Buoyancy-inhibited Ekman flows create a persis- 
tent angular velocity differential between the crust 
and core with Afl/fl as high as 10 _1 0j- 

Once a detection is made, the analysis in section llLTl will 
need to be generalized to distinguish between, and quan- 
tify the relative contribution of, different neutron star 
populations as well as other (e.g. cosmological) emission 
mechanisms. Parameter estimation in this case is signif- 
icantly more complicated and needs to be evaluated in 
the context of specific detection algorithms. We do not 
attempt it here in view of the uncertainties outlined in 
the previous paragraph. Nevertheless, by way of illus- 
tration, a cross-correlation search requires detection over 
a finite frequency band in order to determine the form 
of the power-law. This can be complicated by the pres- 
ence of multiple sub-populations of pulsars, e.g., a signal 
from both accreting and non- accreting systems, or the 
overlap of a signal from another stochastic source. Ad- 
ditional complications arise by noting that the overlap 
reduction function is typically small at low frequencies 
[59l ] . although Advanced LIGO does have some narrow- 
band capabilities. The subject of parameter estimation 
in anticipation of future detectinos is the subject of on- 
going work. 



Appendix A: Stratified Turbulence 

Shear-driven turbulence in a fluid that is stably strat- 
ified against thermal convection is a subtle phenomenon. 
Many open questions persist regarding terrestrial experi- 
ments with Navier-Stokes fluids, let alone exotic superflu- 
ids in neutron stars. A proper treatment of stratification 
lies well outside the scope of this paper. In this appendix, 
we flag some of the key issues briefly and point the reader 
to some useful references, in anticipation of further stud- 
ies. The issues are also canvassed in sections 2 and 3.3 
in Ref. @. 

The latest results on stratified turbulence come from 
large-scale (e.g. 1024 x 512 x 512) direct numerical sim- 
ulations in three dimensions, e.g., [oOl l6l| and references 
therein. The simulations are controlled by two vari- 
ables: the activity parameter, / = e/vN 2 , where N is 
the Brunt- Vaisala, frequency, and the Richardson num- 
ber, Ri = A^ 2 {dvrj,/ dr)~ 2 (or equivalently the reciprocal 
of the squared Froude number). As a rule of thumb, 
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stratified turbulence is three-dimensional for / > 7 [62j ]. 
with a fully developed inertial range and dissipation at 
small scales. It fossilises or "collapses" down to two di- 
mensions for / < 7, whereupon viscous shearing dominat- 
ing vertical momentum transport and dissipation occurs 
at large scales (i.e. there is no inertial range). Figure 
18 in Brethouwer et al. [6(| delineates these two regimes 
(as well as unstratificd turbulence) more precisely in the 
/-Ri ' plane. Given the scalings 



/ = 4 x KV 



Ri = 3 x 10 5 



AO 



v N- 1 / N 



lrads -1 / Vln^s- 1 / V500rads _1 

(Al) 



An \ 2 / n 



lrads -1 / \500rads~ 1 



(A2) 



most neutron stars lie in the three-dimensional regime, 
near the boundary between strong and weak stratifi- 
cation. Buoyancy suppresses radial motions above the 
Ozmidov scale 61 1 



(e/N 



3\l/2 



An 



1 rad s 



3/2 



N 



500 rads" 1 



-3/2 



m (A3) 



Three other scales in the problem are the Corrsin scale, 
e 1 ^ 2 (dv r j ) /dr)~ 3 / 2 , below which anisotropic shear produc- 
tion is weak, the Kolmogorov scale for viscous dissipa- 
tion, (t^/e) 1 / 4 , and the radius of the star. 

Even when the turbulence fossilises for I < 7, the 
flow remains turbulent in spherical shells rather than 
in three dimensions. This phenomenon is observed in 
high-resolution numerical simulations 6^, [6l[ and also 
in the Earth's atmosphere and oceans 63[. Large-scale, 
vertically sheared (i.e. streamwise elongated), horizon- 
tal motions persist under strong stratification, e.g. pan- 
cake vortices in spherical shells, and interleaved laminar 
and turbulent "lasagne" layers. The structures are all 
non-axisymmetric in general. Furthermore, they are in- 
termittent, recurring erratically in bursts even when the 
system is statistically stationary (achieved in the sim- 
ulations by throttling the mean shear (6l1|). Thus, a 
stochastic gravitational wave signal is expected under 
a wide range of stratified conditions, even though its 
character changes with /. Likewise, momentum is trans- 
ported vertically, whether the turbulence is two- or three- 
dimensional. Chung and Matheou [6l[ showed that the 
vertical diapycnal diffusivity scales with the activity pa- 
rameter app roximately as vl 1 ^ 3 divided by the Prandtl 
number [64| . with the details depending on lo/R*- 

How might the predictions in this paper change in light 
of the above stratification physics? The effects enter in 
three places. First, the power spectral density P(k) cx k a 
arguably changes from a = —11/3 (/ > 7; forward cas- 
cade) to a = —3 (/ < 7; reverse cascade). However, C(0) 
changes by less than 20% between these two cases [9[ , and 
the scalings of C(0) and t c with Ail, R+ and M* are the 



same for all a < —5/3 (i.e. the signal is dominated by 
k s ). Second, for I < 7, the signal becomes intermittent 
(see above) albeit still statistically continuous. Again, 
though, this effect is likely to wash out when observing 
a cosmological background in which multiple sources are 
superposed. Third, stratification may invalidate our as- 
sumption of isotropy, e.g., when pancake vortices form 
for / < 7. This is a genuinely open (and difficult) ques- 
tion which deserves further study. It is complicated by 
the fact that even unstratified high-Re turbulence dis- 
plays long-lived coherent anisotropic structures like hair- 
pin vortices, e.g., in terrestrial wind-tunnel experiments 
[65| ; see also Ref . Moreover, the quantum mechanical 
turbulence ("vortex tangle") in a superfluid is polarised 
on large scales even though it is isotropic locally [TH, |66[ 
and anisotropic on intermediate scales due to patchy mu- 
tual friction [12| ■ Ultimately, the magnetic field must also 
be treated, raising other difficult issues. 



Appendix B: Wiener-Khintchine theorem 

The derivation of equation (fT5"]) from equation (TT3]) re- 
quired careful application of the Wiener-Khintchine the- 
orem to ensure that the energy flux does not diverge but 
rather is proportional to the emitting time. Our deriva- 
tion follows closely that outlined by Pottier [Zol ]. 

Consider a stochastic process, X(t), which is both real 
and stationary. Realisations of this process, x(t), are not 
square-integrable [i.e., dt \x{t)\ 2 diverges], as x(t) 
docs not vanish as t — > oo. Hence, we consider a finite 
time interval, T, and define a truncated time scries 



X T (t) 



X(t) 0<t<T 
elsewhere 



(Bl) 



The Fourier transform of Xt is generally defined in terms 
of an integral from — oo to oo and reduces here to 



X T {v) = [ dtX T {t)c 
Jo 



2-nivt 



(B2) 



for the truncated time series. Moreover, the Fourier co- 
efficients are expressed as 

A n = ^ J T X{ty v ^dt - ^X{v n ), (B3) 

where the last relation assumes a fixed T . 

The power spectral density, S{v), is proportional to 
the mean square of the Fourier transform 



S{v) = ±(\X T ( V )\ 2 ). 



(B4) 



Assuming S{v) to be a continuous function of i>, and 
taking the limit T — > oo, one can show 

{X{v)X{v'Y) =2tt5{v -v')S{v). (B5) 
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Evaluating this at v = v 1 leads to an infinite energy per 
unit area per unit frequency: a turbulent neutron star 
emitting for an infinite length of time radiates an infinite 
amount of energy. For a finite emitting lifetime, T, we 
have from (IB3I) 



<KI 2 ) = ^f dt Io dt ' (I(i)1(f ' )>e '"" (i "''' 

(B6) 

The autocorrelation function in the integrand is only a 
function of r = t—t'. Integrating over the square domain 
< t < T, < t' < T, we find 



<KI 2 > 



- 1 /' dr (l |T| » 



and hence from (IB4I) 



T 



(X(t)X(t'))e w ^, (B7) 



S(u) = J t dr \ l U j (X(t)X(t')) c-. (B8) 

Taking the limit as T — > oo yields the standard Wiener- 
Khintchine theorem that S(v) is the inverse Fourier 
transform of the autocorrelation function. 



In the present application, we identify X(t) with 
dhj^/dt. Hence (X(t)X(t')) falls away with r on the 
turbulence decorrelation timescale, r c , which is of the or- 
der of milliseconds, significantly longer than the relevant 
emitting time (~ 10 9 yr). We therefore have T ^> \r\ and 
hence 



(\X(v) 2 \)=T dr {X{t)X{t'))c WT . (B9) 
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